Superscattering emerging from the physics of bound states in the continuum

We study the Mie-like scattering from an open subwavelength resonator made of a high-index dielectric material, when its parameters are tuned to the regime of interfering resonances. We uncover a novel mechanism of superscattering, closely linked to strong coupling of the resonant modes and described by the physics of bound states in the continuum (BICs). We demonstrate that the enhanced scattering occurs due to constructive interference described by the Friedrich-Wintgen mechanism of interfering resonances, allowing to push the scattering cross section of a multipole resonance beyond the currently established limit. We develop a general non-Hermitian model to describe interfering resonances of the quasi-normal modes, and study subwavelength dielectric nonspherical resonators exhibiting avoided crossing resonances associated with quasi-BIC states. We confirm our theoretical findings by a scattering experiment conducted in the microwave frequency range. Our results reveal a new strategy to boost scattering from non-Hermitian systems, suggesting important implications for metadevices.


S1. QNM PERTURBATION THEORY
In this section, we introduce a general perturbation scheme to understand the evolution of QNMs of an arbitrary geometry. Later, we present simplified expressions for the coupling between two QNMs. We utilize the latter to analyze the deformation of a sphere to an ellipsoid of revolution in the vicinity of two QNMs radiating as two different multipoles.
In a coupling experiment, one or more parameters of the cavity are varied. In the perturbed system (p), the QNMs are linear combinations of the QNMs of the unperturbed system (0), i.e.Ẽ (p) µ . Utilizing first order perturbation theory [1] yields a generalized eigenvalue problem for the 'mixing' coefficients C µ , (organized in a ket vector for compactness):Ω (0) |C =ω (p) (U + V (p) )|C , (S.1) whereΩ (0) is a diagonal matrix containing the eigenfrequencies of the unperturbed cavity, and U is the identity matrix. The ± in the superscripts indicate the field is evaluated at the outer (+) or inner (-) boundary of the cavity.
In what follows, they are omitted for brevity. V (p) encodes the effect of the perturbation, and each component is determined by overlap integrals of the unperturbed modes over the surface of the cavity ∂V . δh(r) is the geometrical deformation, where δh(r)n corresponds to a perpendicular shift of the resonator boundary from the undeformed to the deformed surface [1], as illustrated in the insets of Fig. S1(a).
Equations (S.1) and (S.2) are general and can be applied to any open dielectric cavity. This naively seems to imply that arbitrary modal overlaps can occur. However, in canonical geometries, there exist strict symmetry rules that prevent energy exchange between QNMs with different multipolar origin. To illustrate this, consider a sphere of radius R supporting a QNM radiating to the ED channel (a), withẼ (0) a ∝Ẽ p (r), and a low radiative QNM which only couples to the magnetic quadrupole, labeled b, withẼ (0) b ∝Ẽ M (r), (the tilde denotes that the multipole fields radiated by each mode must be evaluated at its complex eigenfrequency). The latter effectively corresponds to a quasi-BIC [2]. We first consider the case of a radial deformation [left inset in Fig. S1(a)], which preserves the spherical symmetry. For a bimodal system, the eigenvalue problem in Eq. (S.1) reduces to The dependence of all V ij elements with the boundary shift has been omitted for clarity. In spherical coordinates, the boundary shift δh(r) = δR, where δR is a constant determining the increase (or decrease) of the sphere radius. Evaluating the perturbation term V ab then reduces to integratingẼ p (r) ·Ẽ M (r) over the whole solid angle. Therefore, due to the orthogonality of the multipole fields, V ab = V ba = 0 in this first case. Then, the eigenvectors are |C u ⟩ = (1, 0) T and |C d ⟩ = (0, 1) T , where T denotes transpose. This effectively means that the fields of the perturbed system will be identical to the unperturbed one (not necessarily the eigenfrequencies, which can be shifted since the V ii terms are nonzero). While we have worked on the specific example of ED and MQ modes, the result is general: since multipole fields are orthogonal on the unit sphere, any perturbation that preserves the spherical symmetry cannot mix modes matched to different multipole channels.
We now consider the case of interest for this work, namely, the situation when the spherical symmetry is broken. We investigate the transition from a sphere to an ellipsoid of revolution [right inset in Fig. S1(a)]. The boundary shift takes a more complicated form, but it can be shown that, if the ratio between the in-plane and out-of-plane axis of the ellipsoid η is close to unity, to second order in η, δh(r) = δh(θ) ≈ R(1 − η) cos 2 θ, where θ is the polar angle in spherical coordinates. The V ab term is: The integral in Eq. (S.4) can be nonzero because (i) the scalar product between two multipole fields of different order does not vanish, and (ii) the function cos 2 θ prevents the direct application of the orthogonality between the multipole fields in the unit sphere. Since the boundary shift does not depend on the azimuthal angle, the structure retains the rotational symmetry in the x-y plane, so that its point symmetry changes from O(3) to D ∞h (in Schönflies notation). Then, it can be shown through group theoretical arguments that V ab does not vanish, since it couples modes radiating u expressed in the basis of the original 'pure' multipolar QNMs of the sphere a, b. When r ⊥ /r ∥ ̸ = 1, V ab ̸ = 0, and therefore the mixing angle is nonzero, inducing hybridization between a, b. (d) Intrinsic ED moment ofẼ which radiates only as a magnetic quadrupole (a quasi-BIC). Mode hybridization leads to leakage to the ED channel once r ⊥ /r ∥ ̸ = 1. Conversely,Ẽ (p) d always retains its dipolar nature, since when r ⊥ /r ∥ = 1,Ẽ as electric and magnetic multipoles of different order [3], which now belong to the same irreducible representation.
Solving Eq. (S.3) we obtain the perturbed eigenfrequenciesω u,d = E ∓ ∆ sec(2ϕ), where: and Similarly, we define the new eigenvectors as: Equations (S.9) and (S.10) confirm that the new QNMs, once the spherical symmetry is broken, are linear combinations of the QNMs in the sphere. In analogy with a two-level model in quantum physics [4], we call ϕ the mixing angle, since it determines the energy exchange between the original QNMs of the sphere. Making the replacement r = ω b /ω a , Eqs. (S.9) and (S.10) correspond to the hybrid resonances (matched to more than one multipole) discussed in the main text. From the expression of the perturbed eigenfrequencies, we also clearly observe the formation of an avoided crossing, which ultimately depends on the coupling strength Notably, if V ab = 0, then ϕ = 0 and there is no mixing, since again we recover |C u ⟩ = (1, 0) T and |C d ⟩ = (0, 1) T . However, in all other situations there is hybridization. This indeed corresponds to the Friedrich-Wintgen mechanism [5].
We are now able to delve into the physics of a realistic nanocavity. In particular, we study the axial deformation of a Si nanosphere with radius 100 nm, discussed in the main text . In Figs. S1(a) and S1(b) we have plotted the real and imaginary parts of the eigenfrequencies recovered with the two-level model in Eq. (S.3). The deformation is parameterized by the ratio between the out-of-plane and in-plane axis of the ellipsoid, r ⊥ /r ∥ = 1/η.
As expected, we observe an avoided crossing of the real part, as well as a reduced imaginary part for mode u when r ⊥ /r ∥ = 1, corresponding to the quasi-BIC condition. In accordance with the two-level picture described in the previous paragraphs, Fig. S1(c) shows how the eigenvector |C u , once the spherical symmetry is broken, receives contributions from the two original QNMs. Most importantly, due to the hybridization, mode u can now radiate as an ED. This can be seen in Fig. S1(d), where we have plotted the ED moment of each QNM (evaluated at their eigenfrequency). For QNM u, the ED moment is zero for a sphere, due to the quasi-BIC condition, while QNM d has a large ED moment. Breaking the spherical symmetry inmediately results in coupling between the QNMs: QNM u, due to the hybridization with d, now possesses an ED moment. This mechanism is essential to overcome the single channel limit, as explained in the main text.

S2. CONNECTION BETWEEN PERTURBATION THEORY AND TCMT
In this section we outline the link between QNM perturbation theory and TCMT, and the approximations implied. Based on the link, we obtain suitable analytical forms for the evolution of the TCMT parameters discussed in the main text. A connection can be made by writing the equation for one perturbed QNM, from Eq.(S.3), and comparing it with its counterpart TCMT equation:ω In the previous, we have denoted byω the eigenfrequency of the new QNM resulting from the coupling between |a, b⟩. The superscript (1) used in the TCMT equation will become clear in what follows. Dividing the first equation by 1 + V aa leads toω a /(1 + V aa ) represents the evolution that QNM |a⟩ would have in the absence of coupling with |b⟩. Comparing the TCMT equation with the above, it appears that the coupling coefficient can be defined as The dependence onω can be neglected under the assumption of weak coupling inherent to TCMT. Namely, the effect of κ ij (ω) can only be felt in the vicinity of the degeneracy, whereω b )/2 and δ is a small number, which can be neglected, yielding Thus, we can now build the effective Hamiltonian of TCMT as (S. 16) We now impose the second main approximation of TCMT and assume high-Q modes, so that κ ab = −κ * ba . Further, we assume all coupling takes place in the near field, so that κ ab = κ ba = −iκ ∈ iR. With this, H 0 takes the simple form given in Eq.(2) in the main text: (S.17) Finally, to study the Friedrich-Wintgen mechanism, based on the arguments above, we consider the following functional forms:ω where the a i are constants, representing the surface integrals within the V ij , and ζ is the perturbation, extracted from δh(n), assumed to be independent on the spatial coordinates.

S3. DERIVATION OF EQ.(2) OF THE MAIN TEXT
For a given R-matrix, it can be shown that the total scattering cross section stems from the contributions of all scattered multipoles, Every σ τ can be found as [6,7]: where R τ ν is an element of the reflection matrix R(ω) defined in the main text, δ τ ν is the Kronecker delta and |a ν | 2 is the power carried by the plane wave in the multipole ν [8]. |a ν | 2 /I 0 = (2ℓ ν + 1)λ 2 /2π, where ℓ ν is the angular momentum associated with the multipole ν. We consider a hypothetical scenario where channel 1 (for instance the electric dipole), receives contributions from another channel τ = 2. For clarity, we normalize Eq.S.21 by |a 1 | 2 = 3λ 2 /2π = σ Max (the limit to the channel cross section in the absence of coupling). The scattering cross section for channel 1 is then given by: where we have assumed, for simplicity in the toy model, that the power carried by the plane wave in both channels is the same (i.e. a 1 = a 2 ).

S4. DERIVATION OF f (ω)
In this section, we derive expressions for f (ω) for the two cases considered, (i) single mode approximation, (ii) two mode approximation. We start with (i). Using TCMT, we write the following equations describing the system dynamics: In Eqs.S.24-S.25, the d j are the coupling coefficients of the resonant QNM with the incoming waves j. They are related to the radiative losses on each channel as γ j = d 2 j , and the total radiative loss is given by γ = j γ j . With a time dependence in the form e −iωt , the R-matrix can be found from with M = dd T , and d is a column vector containing the d j coefficients. The off-diagonal elements R ij are nonzero due to the term M 12 = d 1 d 2 = √ γ 1 γ 2 of the matrix M. The function f (ω) in Eq.S.23 can be evaluated as: Next, we consider the system is driven by two resonant QNMs, which corresponds to case (ii). Each one is matched only to a given channel τ = 1, 2. However, the two modes can also interact with each other in the near field, with a coupling constant denoted by κ. Altogether, the TCMT equations can now be formulated as: Eqs.S.28-S.31 can be written in a more compact form: where A is a column vector containing (|a⟩ , |b⟩) T , s ± is a column vector containing the incoming (or outgoing) amplitudes, H 0 is an effective Hamiltonian governing the mode evolution, and D is a diagonal matrix containing the coupling coefficients of the modes to the scattering channels, which reads: The d j are connected to the decay rates of each QNM as γ a = d 2 1 and γ b = d 2 2 . If the structure undergoes a suitable geometrical perturbation parameterized by ζ, the QNMs can couple. Within the framework of first order perturbation theory, the effective Hamiltonian H 0 in the basis of the original QNMs reads: the coupling constant κ ∈ R, as discussed in section S2. The new QNMs of H 0 can be expressed in the basis of the unperturbed ones as |u, d⟩ = (∆ ± G) |a⟩ + κ |b⟩, where 2∆ =ω a −ω b and G 2 = E 2 + κ 2 , with 2E =ω a +ω b . The new fields are necessarily a linear combination of the unperturbed ones. However, the 'hybrid' QNMs no longer diagonalize the D matrix: each QNM is matched to the two channels, in a proportion dictated by κ. Hereon, it is assumed that γ b << γ a . Moreover, the channels and the modes are in the basis such that radiative losses are minimal for |b⟩ at κ = 0. Mode |b⟩ is therefore a quasi-BIC.
Therefore, in the presence of coupling between two QNMs matched to different multipoles, the R-matrix is no longer diagonal. Following the same rationale that led to Eq.S.27, we reach the final expression for f (ω): . (S.36)

S5. DEPENDENCE OF THE SCATTERING CROSS SECTION WITH SMALL OHMIC LOSSES
As discussed in the main text, the ability to control the Q-factor of the superscattering resonance is a key feature of a super multipole, since it allows to maintain a large scattering cross section even in the presence of ohmic losses. Indeed, large Q-factors come at the disadvantage of stronger dissipation, leading to a steeper drop in the scattering cross section. To demonstrate this, consider a resonance matched to a single multipole τ = 1, having the normalized scattering cross section (at resonance): In the previous, γ 0 denotes the loss rate induced in the mode by absorption. For small absorption, γ 0 can be taken to be proportional to the imaginary part of the refractive index 'δ'. We estimate the rate of change of the scattering cross section with increasing δ as: where Q = ω 1 /2γ 1 is the original Q-factor in the absence of ohmic losses, and ω 1 is the resonance frequency. A similar expression can be derived in the presence of leakage to multiple radiation channels.

S6. CARTESIAN EXPRESSIONS OF THE MULTIPOLE COEFFICIENTS
It is sometimes more convenient to express each scattering channel τ in the Cartesian basis [9]. We emphasize that each multipole moment (whether expressed in the spherical or the Cartesian basis) constitutes a scattering channel, and the contribution to the cross section in a sphere is limited to the usual bound. Recently, elegant formulae for the multipole moments in the Cartesian basis in terms of the source currents have been proposed [9]. We provide them here for consistency, with a modification of the prefactors to account for power normalization: The prefactors in Eqs. (S.39)-(S.42) are: The integrals in Eqs. (S.39)-(S.42) are performed over the resonator volume. The polarization currents P are defined as P = δε(r)E(r), where δε(r) = ε(r) − ε 0 is the permittivity contrast between the resonator and free space, and E(r) is the total field within the resonator. For clarity, the spatial dependence of P has been omitted everywhere.

S7. QNM EXPANSION OF THE MULTIPOLE COEFFICIENTS
QNMs form a complete basis for the total field inside the resonator [10]. Consequently, we can write: where θ(r) is a function that equals one within the volume of the resonator V , and zero outside. There are several equivalent forms for the expansion coefficients α µ (ω) [10,11]. Assuming a non-dispersive material, normalizing the QNMs with the Sauvan norm [12], and making use of reciprocity arguments, the α µ (ω) are given by [10]: The background field is, in our case, a normally incident plane wave, e.g. E b (r; ω) = E 0 e ikz e x . We can now define the QNM contributions to the total polarization current asP µ = δεẼ µ , so that P = µ α µPµ . Substituting this expression into the multipole moments in Eqs. (S.39)-(S.42) yields a rigorous QNM expansion of the multipole moments at real or complex frequencies.

S8. EXPERIMENTAL METHODS
A Taizhou Wangling TP-series microwave ceramic composite is used as a dielectric material for resonators (disks) fabrication. According to the manufacturer data sheet, relative permittivity of this ceramic is 22 ± 1 and loss tangent is tan δ ≈ 1 · 10 −3 at 10 GHz. A set of disks is fabricated from the ceramic plates with the use of a precise milling machine. The radius of disks is 4 mm. The sample to be measured is composed of several disks to obtain several resonators with thicknesses of 6.0, 6.9, and 10.0 mm. The sample is mounted on a support made of ROHACELL ® 71 HF plate, whose relative permittivity is 1.01 ± 0.05.
To measure the total extinction cross section, the sample is placed in an echoic chamber and illuminated by a normally incident linearly-polarized waves radiated and received by a pair of HengDa Microwave HD-10180DRA10 horn antennas. A photograph of the experimental setup is shown in Fig. S2(a). The antennas operating range is 1-18 GHz. The distance between each horn antenna and the sample is fixed at 2.0 m. The antennas are connected to the ports of Rhode & Schwarz ZVA50 Vector Network Analyzer (VNA) by 50-Ω coaxial cables. Using VNA the S21-parameter (transmission coefficient) is detected and analyzed by a special computer software. Forward scattering is obtained from the measured transmission coefficient. The total extinction cross section is extracted from the measured complex magnitude of the forward-scattered signal by means of the optical theorem [13]. A stainless steel sphere with the radius 40.00 ± 0.04 mm and the corresponding analytical Mie solution are used for the signal calibration. To reduce unwanted reverberations between antennas the time gating technique is used. Then a standard filtering procedure is applied to the raw data to remove the background noise.
A photograph of the measurement setup for mapping the electric near-field is presented in Fig. S2(b). The transmitting antenna generates a linearly polarized quasi-plane-wave whose polarization is orthogonal to the symmetry axis (the z axis) of the disk resonator. The receiving antenna is substituted by an electrically small probe. This probe is a Hertz dipole which is oriented along the polarization direction of the incident wave. The probe detects the dominant component (E x ) of the scattered electric near-field. A LINBOU near-field imaging system is used for the near-field mapping. The scan area around the sample has dimensions 50 × 50 mm 2 in the x-z plane. In the measurements, the probe automatically moves in the scan area with a 1-mm step along two orthogonal directions. At each probe position, both the amplitude and the phase of the scattered electric field are measured. The exception is the area near the sample. This area is experimentally inaccessible because of the technical limitation on the distance between the moving probe and the vertically standing sample. In the next section we provide further details on the numerical simulations of the ceramic rods utilized in the experiment.  Here, the electric dipole can exceed the single channel limit, and there is no need for crossing of the resonant frequencies.

S11. ED RESONANCES IN THE SPHERE AND THE ELLIPSOID
In this section we compare the scattering cross section of the first ED resonances in a Si sphere (n = 4) vs. an ellipsoid, to further demonstrate the correctness of the single-channel limit in the main text (Fig.S5). For the sphere, all the resonances are always bounded to the single channel limit. However, the optimized ellipsoid can greatly surpass this limit. Moreover, the maximum enhancement can be clearly seen to occur at the super ED resonance described in the main text.
FIG. S5. Lowest order ED resonances of a Si sphere and a Si ellipsoid (refractive index fixed at n = 4), with aspect ratio b/a = 0.8 vs. a/λ. b, a are the semiaxis of the ellipsoid (or the radius of the sphere). a is kept fixed at 100 nm. For the sphere, all the resonances are bounded to the single channel limit. However, the optimized ellipsoid can greatly surpass this limit. Moreover, the maximum enhancement can be clearly seen to occur at the super ED resonance described in the main text. Only the ED scattering cross section is plotted.

S12. COMPARISON WITH PREVIOUS WORKS ON SUPERSCATTERING WITH NON-SPHERICAL SCATTERERS
In this section we study in more detail previous works [14,15] investigating non-spherical scatterers with enhanced scattering, and demonstrate that they do not display super multipoles. First, we consider the cylindrical resonator investigated in Ref. [14]. In the latter, the authors presented a Ag nanocylinder embedded in a polymer with refractive index 1.4, and claimed to reach the superscattering regime. We note that all the cross sections in their work are normalized to arbitrary units, making it difficult to determine the scattering enhancement. In consequence, we have reproduced their simulations with the help of COMSOL multiphysics. Fig.S6 shows our calculations of the scattering cross section using the same parameters (details are given in the caption). The calculated scattering cross section is only slightly above the single-channel limit (approximately 1.2 times). We attribute the small efficiency to the presence of Ohmic losses from Ag. Most importantly, the mild enhancement is due to the spectral overlap of the ED and EQ resonances, which are associated to two distinct eigenmodes of the cylinder [ Fig. S6(b)]. Each multipole, by itself, does not exceed the single-channel limit. Hence, despite investigating a similar geometry, the authors did not observe the same effect that we are describing in our work, where a single multipole exceeds the limit.
Next, we consider the structures investigated in Ref. [15]. The authors considered Ag ellipsoids of extreme subwavelength sizes, in the order of 6% the wavelength, [as depicted in the inset of Figure S6(c)]. We have calculated the scattering cross section of one of such ellipsoids near their localized plasmon resonance, [ Figure S6(c)]. The scattering cross section is 1000 times below the single channel limit. Therefore, these ellipsoids are very weak scatterers. Moreover, Figure S6(d) compares the magnitudes of the scattering and absorption efficiencies. Absorption is almost 20 times larger. To understand these results, we should clarify that the authors in Ref. [15] were not interested in maximizing the scattering cross section of a single particle; their figure of merit was the extinction per unit volume. In addition, as we have now verified with full-wave simulations, the response of such particles will always be strongly dominated by absorption, in contrast to our dielectric nanoparticles.  [14]. Following the latter, the nanoparticle is described by the Drude model ε(ω) = ε∞ − ω 2 p ω(ω+iγc) , with the following parameters: ε∞ = 4, ωp = 1.336 · 10 16 rad/s, γc = 1.108 · 10 14 rad/s . The scatterer is embedded in a host environment with refractive index n h = 1.4 , hence the single channel limit is given by σMax = 3 3πn h λ 2 . The diameter of the cylinder is 80 nm, and the height is 120 nm. All corners have been rounded to avoid unphysical field singularities. (a): Multipole decomposition of the scattering cross section. Inset: scheme depicting the illumination setup. The shaded green area indicates the spectral region with maximum scattering enhancement. The enhanced scattering occurs due to the overlap of the ED and EQ contributions, following the conventional mechanism of Fan et al [6,8]. Note that none of the two multipoles exceeds the single channel limit, and the total cross section is only slightly above it. (b) Upper panels: field distributions of the two involved QNMs. The rightmost panel is the ED QNM, and the leftmost is the EQ QNM. Lower panels: Evolution of the surface charge during half an oscillation period at the wavelength with maximum scattering enhancement. (c)-(d) Ag ellipsoids investigated in Ref. [15]. (c) Multipole decomposition of the scattering cross section of an ellipsoid composed of Ag (material data obtained from Johnson and Christy [16]), with semi-major axis 10 nm. The scattering cross section is entirely dominated by the ED, and shows values almost 1000 times below the single channel limit. Inset: scheme of the illumination setup and the scatterer under consideration. (d) Comparison between scattering and absorption efficiencies (in arbitrary units). The absorption cross section is almost 20 times larger than the scattering cross section.